Coulomb blockade in electron transport through a Ceo molecule from first principles 
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We present results of spin-unrestricted first-principles quantum transport for a gated Ceo molecule 
weakly contacted to Al electrodes, making emphasis on the role played by the electronic localiza- 
tion and the spin degree of freedom. As expected, the conductance presents a series of peaks as 
a function of a gate voltage, demonstrating that transport in the Coulomb blockade regime can 
be properly treated within a first-principles scheme. A well-known manifestation of the interplay 
between Coulomb interaction and the spin degree of freedom in atoms and molecules, the Hund 's 
rule, determines the sequence of conductance peaks. 
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I. INTRODUCTION 

Since molecular electronical is the focus of attention on 
the part of investigators in search of the ultimate device 
miniaturization, many attempts have been made to con- 
tact molecules to metallic electrodes. With the appropri- 
ate end groups these molecules can bind to the electrodes 
and form fairly stable molecular bridges as initially shown 
by Reed et al^. Since then, a combination of chemi- 
cal functionalization and nanofabrication techniques are 
being pursued with different degrees of success^. Simi- 
larly, the road to describe and understand theoretically 
how these contacted molecules carry current is plagued 
with technical and conceptual difficulties^. One of these 
difficulties stems from the use of density functional the- 
ory (DFT), on which most of the codes developed for 
the calculation of "first-principles" transport are cur- 
rently basec&£iL2i2iiii. The notorious discrepancies be- 
tween experiments and theory are being attributed by 
some groups to the failure of standard DFT in describ- 
ing the behaviour of the quasiparticles at the Fermi level 
when charge density variations are pronouncecUiiSii^. 

While failures attributed to the use of DFT might be 
critical on cases where strong charge localization occurs 
and the basis for this criticism is worth pursuing further, 
it is also imperative to understand the role played by all 
the others factors at play. One of the factors that has not 
been fully (not even partially) explored is the spin, ft is 
quite obvious that magnetism is not a lesser factor when 
it comes to device functionality. While a large amount 
of experimental work has been done in large area mag- 
netic junctions and magnetic multilayers due to the huge 
technological impact of these systems^, a deep theoreti- 
cal understanding of electron transport in these systems 
is still lacking. Magnetic nanocontacts, for instance, ex- 
hibit a very rich and complex behavior which, at this 
moment, is far from being understoocU£*i£ii!. Some non- 
magnetic metallic nanocontacts are also expected to ex- 
hibit weak magnetism due to the low coordination of the 
atoms at the breaking points. Magnetic molecules are 
also the subject of several experimental studies from the 
transport point of viewi^. Furthermore, even if all the 






FIG. 1: (a) A molecular bridge formed by a Ceo molecule 
weakly attached to Al electrodes with pyramidal form. 
Coulomb blockade is expected to occur due to the formation 
of quasi-localized molecular orbitals. (b) An Al nanocontact 
where electronic localization is not expected. 



elements composing a molecular bridge, namely, the elec- 
trodes and the molecule, are not intrinsically magnetic, 
charge density localization due to poor coupling of the 
molecule to the electrodes can help localize a spin den- 
sity at the molecule. This last scenario is our focus in 
this work. 

To date, most first-principles codes for electronic trans- 
port do not take into account the spin degree of free- 
dom since this adds an important degree of difficulty to 
the well-known computational comple xity of thi s type of 
calculations (for exceptions see Refs. |!0l20l2lj ). Here, 
we present results obtained from the spin-unrestricted 
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version of our GAUSSIAN03-based2£ ab initio code for 
quantum transport named ALACANT (ALicante Ab ini- 
tio Computation Applied to NanoTransport) 23 . For illus- 
tration purposes, we consider here a prototypical molec- 
ular bridge: A Cqo molecule (weakly) contacted to Al 
electrodes and in the presence of a nearby metallic gate 
[see Fig. ^a)]. Transport through a Cgo molecule has 
been previously studied theoretically for Al and Au elec- 
trodes with first-principles methods, but without consid- 
ering separate spin transport channels^* 8 * 2 ^ 2 ^. We ex- 
plore here the importance of the spin channel separation 
as well as that of recalculating the conductance as the 
Fermi energy is changed by a gate voltage V g . If the Ceo 
molecule is weakly contacted to the electrodes, the charge 
at the molecule gets strongly localized and transport oc- 
curs in the Coulomb blockade regime. Despite of bulk Al 
being nonmagnetic, spin is relevant for transport through 
the weakly bonded Ceo molecule. The ground state elec- 
tronic configuration of an isolated Cgo molecule is non- 
magnetic, with a three-fold degenerate lowest unoccupied 
molecular orbital (LUMO), t u , and a five-fold degenerate 
highest occupied molecular orbital (HOMO), h u . Charge 
transfer from the electrodes or changes in the Fermi en- 
ergy induced by nearby gates can change the charge state 
of the molecule and, therefore, the molecular electronic 
configuration can become magnetic. In general terms, 
transport through weakly contacted molecules presents 
zero-bias conductance peaks as a function of V^^SLS^. 
The molecule gets filled or emptied one electron at a time 
by the action of this voltage, opening and closing shells 
(almost) succesively. Thereby the necessity of consider- 
ing the spin degree of freedom in these situations. In 
addition to that, we find that Hund 's rule determines 
the sequence of conductance peaks, something that has 
been observed in poorly contacted nanotubes by several 
group a 2 ?' 28 . To our knowledge, this is the first time that 
transport in the Coulomb blockade regime is fully ana- 
lyzed from first principles. 



II. THEORETICAL BACKGROUND 

A. Green's function formalism for open-shell 

systems 

The design and fabrication of electronic devices at the 
molecular and atomic scale has posed new challenges to 
theorists. The basics to calculate the zero-bias, zero- 
temperature conductance, G, in a molecular bridge or 
metallic nanocontact were established by Landauer long 
before the concept of nanoelectronics was commonplace. 
In Landauer 's formalism Q is proportional to the quan- 
tum mechanical transmission probability for the electrons 
at the Fermi energy, Ep, to cross the molecular bridge 2 ^: 

g=^(E F )+T l (E F )}. (1) 



In this expression the contributions from spin up (f) 
and spin down (j) channels have been explicitly sepa- 
rated while the contribution from all the orbital chan- 
nels has been condensed in T. For simplicity we ne- 
glect spin mixing due to spin-orbit scattering and ex- 
clude the possibilty of noncollinear spin densities, i.e., 
S z is a good quantum number. It is well-known that 
the detailed atomic, electronic, and magnetic structure 
at a nanocontac1*i£ is important and, in order to achieve 
a quantitative level of agreement with experiments, one 
has to rely on first-principles or ab initio calculations, 
typically at the DFT level. For molecular bridges like 
the one in Fig. ^a) this necessity becomes imperative 
due to the impossibility of predicting (i) the broaden- 
ing of these orbitals due to the coupling with the elec- 
trodes and, most importantly, (ii) the positioning of the 
Fermi level with respect to the HOMO and LUMO of 
the molecule. Two different approaches to this problem 
can be found in the literature: One based on calculating 
scattering wave functions^S^i and the other based on 
Green function techniqueo 5 i 8 i 9 i 32 i 33 i 34 i 35 i 36 i 37 , sometimes 
combined with the Keldysh formalismZi^SiiS. The basics 
of our approach, which we extend below to include the 
spin degree of freedom, have been presented in previous 
publications 8 ! 32 ! 38 . 

The main advantage of the numerical implementation 
that we have developed is that relies on a standard quan- 
tum chemistry code such as GAUSSIAN03. This code is a 
versatile tool to perform first-principles or ab initio cal- 
culations of molecules or clusters that incorporates the 
major advancements in the field in terms of functionals, 
basis sets, pseudopotentials, etc.. The way to proceed 
is as follows: A standard self-consistent field electronic 
structure calculation of the molecule that includes a sig- 
nificant part of the electrodes [like that shown in Fig. 
^a)] is performed. This calculation is usually performed 
at a DFT level in any of its multiple approximations. 
The use of configuration interaction or multiple Slater de- 
terminant methods has been recently proposed to study 
regimes of transport beyond the scope of DFT such as 
the Kondo effect^ or simply to improve existing DFT 
results^, but these approaches limit the calculations to 
very simple systems and present some conceptual diffi- 
culties not fully resolved. 

As far as transport is concerned, the self-consistent 
hamiltonian Hjq) (or Fock matrix F^)) of the central 
cluster contains the relevant information. The retarded 
(+) and advanced (-) Green's functions associated with 
the Fock matrices are defined in a standard manner for 
both spin up and spin down species: 

G\f l) =[(E±iS)l-F ni) ]- 1 (2) 

In this expression 1 is the identity matrix and 5 is an 
infinitesimal quantity (in practice we set it to 1CP 10 eV). 
One of the many advantages in the use of Green func- 
tions is that one can incorporate the rest of the infinite 
electrodes in the calculation in a very elegant and simple 



way: 



(El -F m) - 3^(25) 



where 
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J RT(1) 



-(±) 
J LT(i) 



(3) 



(4) 



and Sr (i?)[SL TU) denotes a self-energy matrix 
that represents the spin up(down) hamiltonian of the 
right(left) semi-infinite electrode that has not been ex- 
plicitly included in the calculation. With non-orthogonal 
basis sets, which are the ones implemented in GAUS- 
SIAN03, it has become custommary to use the following 
expression for the Green function: 



(ES - F m) - £ T( n(£0 



(5) 



where the overlap matrix S essentially substitutes 1. The 
self-energy matrices can only be explicitly calculated in 
ideal situations. We describe the bulk electrode with a 
parametrized tight-binding Bethe lattice model with the 
coordination number and effective parameters appropri- 
ate for the type of electrodes. The advantage of a Bethc 
lattice resides in that it reproduces fairly well the bulk 
density of states of most commonly used metallic elec- 
trodes and, at the same time, it does not represent a 
specific atomic arrangement which is also desirable; for 
electrodes are polycrystalline. The expression for the 
Green function in a non-orthogonal basis set presented 
in Eq. [5] is widely accepted in the literature, although it 
does not strictly correspond to the representation of the 
Green function operator in a non-orthogonal se1*22. It is, 
however, convenient since the density matrix P-m) can 
be obtained from it in the standard way: 



TU) 



Im 



dE. 



(6) 



One can now calculate the total electron charge of the 
cluster: 



iV = Tr[(P T +P ; )S]. 



(7) 



In Eq. [7| the trace runs over all the atomic orbitals or 
localized functions composing the basis set. In Eq. [B] 
E-p is determined by fixing the total charge in the clus- 
ter, Q = N + N ion , where Ai on is the fixed ion charge. 
In practice Q needs to be distributed between both spin 
species to comply with the thermodynamic constrain of 
a unique Fermi level for both spin channels. We force 
GAUSSIAN03 to evaluate F T and Fj, using the density 
matrices obtained from Eq. I^and repeat the process until 
convergence is achieved. The spin-dependent transmis- 
sion probabilities that appear in Eq. n can be calculated 
through the well-known expression: 
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FIG. 2: Sequence of conductance curves for different values 
of the total charge Q in the nanocontact shown in Fig. Q^b). 
The curves have been shifted along the energy axis to make 
them coincide as much as possible. 



where, again, the trace runs over all the orbitals and 
Pl,r(E) are twice the imaginary part of the self-energy 
matrices. Finally, in order to single out the contribution 
of individual orbital channels to the current, one can di- 
agonalize the matrix product in Eq. [S] While the size of 
the resulting product matrix in brackets can be as large as 
desired, the number of eigenvalues with a significant con- 
tribution, i.e., the number of conducting channels, will be 
typically much smaller. 



B. Conductance as a function of gate voltage 

Although Eq. [S] gives the transmission as a function 
of energy, only the value of T at the Fermi energy is, in 
general, meaningful (within the limitations inherent to 
DFTJiiiii) . For instance, if a metal plate is placed in the 
vicinity of the molecular bridge and a gate voltage V g is 
applied to it, the total charge of the molecular bridge, 
Q, changes in response to this voltage, as in a classical 
capacitor. Noticing that the gates are typically much 
larger than the molecule, changing the voltage changes 
the net charge of the overall molecule+electrodes system 
and not only of the molecule^. The relation between the 
magnitude of V g and the total charge accumulated in the 
bridge is determined by the capacitance, and this is deter- 
mined, in turn, by the geometry of the whole molecular 
bridge-metallic gate set-up. Having knowledge of this ca- 
pacitance can be of interest, but it is not relevant for our 
purpose here. Therefore, we simply use the net charge of 
the cluster as the independent variable. The whole curve 
T(E) thus needs to be recalculated for each new value of 
the total charge, which, in turn, determines a new Fermi 
energy. In general, T(E) for E ^ E-p does not represent 
a true measurable quantity. Naively, though, one can be 
tempted to describe the effects of this charging as a rigid 
shift of the density of states (DOS) and, consequently, 



T(E) = Tr[r L (£)G«(£)r R (£)G(-)(£0], (8) 



T(E,V g )=T(E-<yeV g ), 



(9) 



where 7 represents a scaling factor associated with the 
specific geometry. This way of effectively accounting for 
the gate voltage is valid for bridges with a very smooth 
DOS such as that in purely metallic nanocontacts as the 
one shown in Fig. ^b). Figure [2] shows T(E) for various 
values of Q. In this case the results essentially confirm 
the naive behavior expected from Eq. (Q- On the other 
hand, if localized or quasi- localized states (e.g., molecular 
states) cross the Fermi energy in performing the energy 
shift, Eq. J3J no longer can be applied. Coulomb repul- 
sion pushes up or down all the empty states, depending 
on whether a molecular orbital fills up or empties. This is 
typically the case when the molecule is loosely attached 
to the electrodes. Charging effects on the molecule dras- 
tically change T(E) as the Fermi energy varies even by 
small amounts. The next section illustrate what, in gen- 
eral, can be stated as 



T(E,V g )^T(E-jeV g ). 



(10) 



III. 



COULOMB BLOCKADE IN TRANSPORT 
THROUGH A C 60 MOLECULE 



Figures and 0] show the zero-temperature transmis- 
sion T as a function of the scattering energy for the 
molecular bridge shown in Fig. Qfa). We have cho- 
sen a pyramidal model for the electrodes in order to 
have a weak contact to the molecule. Similar atomic 
structures for the electrodes close to the contact with 
the molecule are, nevertheless, expected to form in most 
break-junction experiments^. Our present calculations 
are based on DFT theory, but making use of the B3LYP 
hybrid functional 2 ^. The basis sets and pseudopotentials 
are those of Christiansen et al4^£. This combination 
of basis set and functional has been proved to give good 
results at a reasonable computational cost in a variety of 
situations 8 ' 32 - 38 ' 46 . The accuracy of our calculations can 
be systematically improved by resorting to larger basis 
sets, but this is not essential for our discussion here. The 
use of the B3LYP functional provides us, in addition, 
with a semi-phenomenological way of dealing with the 
well-known problems of local density approximations to 
the true functional. The exact non-local Hartree-Fock 
potential is included in this functional, partially taking 
care of the self-interaction problem and inherent dereal- 
ization of charge typical in local or quasi-local approxi- 
mations. This way the value of the transmission off res- 
onance and the shape of the Coulomb blockade peaks is 
expected to approach to the true ones compared to the 
results obtained with LDA or even GGA. 

The six panels in Figs. and 0] present T for both spin 
up (solid lines) and spin down (dashed lines) channels for 
increasing Fermi energy or increasing charge (top to bot- 
tom). It has been shown that a Ceo molecule contacted 
to Al electrodes gets charge transferred from the elec- 
trodes in equilibrium and without the effect of any gate 
voltages 8 . This is a purely chemical effect. In panel (a) 
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FIG. 3: Sequence of transmission curves as a function of the 
overall charge in the whole cluster. The electronic charge in- 
creases from top to bottom, starting from a positively charged 
cluster: (a) Q = 10, (b) Q = 4, and (c) Q = 3. 



of Fig. |3| the Fermi level has been shifted down with re- 
spect to the neutrality point, corresponding to a slightly 
discharged molecular bridge. There the Fermi level lies 
in the HOMO-LUMO gap of the molecule so that the 
bare molecule remains neutral and nonmagnetic. In this 
panel the transmission reveals two degenerate spin chan- 
nels with a series of peaks corresponding to molecular 
states broadened by the coupling to the electrodes. The 
conductance peaks below the Fermi energy reflect the 
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FIG. 4: Same as in Fig. [3] for further increase of the total 
electronic charge: (a) Q = —3, (b) Q = —4, and (c) Q = —10. 



1.0 eV for both spin channels. More molecular states, 
those termed t g , also reflect in the transmission at higher 
energies. 

As the Fermi energy is increased and the doubly- 
degenerate states start to fill up [see Fig. 13(b)], the spin 
degeneracy also begins to be lifted. This broken spin 
symmetry is a consequence of our mean field level de- 
scription. Spin symmetry preserving calculations, which 
may give rise to Kondo physics at low temperatures, are 
out of the scope of our present methodology. It must 
be understood in what follows that the shape and height 
of the peaks would be changed by the Fermi distribu- 
tion function above the Kondo temperature^, which can 
be tuned down to zero by decreasing the coupling of the 
molecule to the electrodes. (We have selected a case with 
a relatively strong coupling for clarity). On further in- 
creasing the net charge [panel (c)] and, consequently, the 
Fermi energy, the spin degeneracy is fully removed and 
the orbital degeneracy of the doubly-degenerate t u or- 
bitals is also lifted as one of them charges up. Again, 
this broken symmetry is a consequence of the mean-field 
approximation. Notice that the spin and orbital sym- 
metry breaking brings the maximum of the conductance 
peaks down to the quantum of conductance e /h, as ex- 
pected for resonant tunneling through an orbital sym- 
metrically coupled to both electrodes (in the absence of 
the Kondo effect). Notice, finally, how the second molec- 
ular orbital that gets filled has the same spin as the first 
one. This is what is expected from Hund 's rule and 
agrees with what various calculations in the literature 
have predicted for Ceo negative ions^. As the Fermi en- 
ergy increases further, the doubly-degenerate spin-down 
orbitals get their degeneracy removed by the charging of 
one of them while the Coulomb interaction pushes the 
other up in energy [see Fig. HJa)]. From here on [Figs. 
0Jb) and (c)] , the remaining spin up orbital and the other 
two empty spin down orbitals get filled one at a time. We 
would like to conclude our analysis by noting that the 
specific order in which the t u orbitals get filled depends 
very much on the bonding geometry to the electrodes and 
the type of electrodes. Also the reliability of the B3LYP 
approximation to DFT in describing the transmission at 
the Fermi energy in the Coulomb blockade regime needs 
to be tested more thouroghly. More sistematic work on 
Coulomb blockade in molecular bridges is obviously re- 
quired. 



DOS corresponding to occupied states and those above 
to empty ones. From now on we focus on the empty ones: 
The three spin-degenerate t u molecular states. For the 
isolated and neutral molecule these states are fully de- 
generate. In panel (a) of Fig. |3ve see that the coupling 
to the electrodes not only broadens these levels, but also 
lifts partially the degeneracy of them. The symmetry of 
these states combined with the chosen bonding geome- 
try to the metal contacts splits these three levels in 2+1. 
This reflects in a transmission peak reaching a value of 
two around 0.6 eV and one reaching a value of one around 



IV. CONCLUSIONS 

We have studied a prototypical example of a molecu- 
lar bridge: a Ceo molecule contacted by Al electrodes. 
It has been shown that spin-unrestricted ab initio trans- 
port calculations in weakly contacted molecular bridges, 
where localization of charge is expected to occur, can re- 
produce the main features of the Columb blockade phe- 
nomenon. In regards to this, it has been shown the ne- 
cessity of recalculating transmission curves as the Fermi 
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energy changes in systems that exhibit charge localiza- 
tion. Transmission curves calculated at different Fermi 
energies can be completely different due to charging ef- 
fects. In addition, Hund 's rule has been shown to play 
an important role in the transport characteristics. 
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